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This paper is being submitted to the Journal of Geophysical Research. It deals with basic 
radiative transfer theory. In particular, it presents a method for using relatively simple 
statistical methods (known as Monte Carlo methods) to perform 3-dimensional radiative 
transfer computations through arbitrarily shaped volumes. This main emphasis of this 
paper is on a technique to deal with non-spherical particles. Such particles, which we 
know to exist in raining clouds, are unique in that they are able to change the polarization 
state of radiation. Because of this, much more sophisticated computations than normally 
used are needed to predict the amount of radiation reaching a satellite. While such methods 
have been described in the literature, this work is unique in that the Monte Carlo method 
really allows the user to arbitrarily specify the size and shape of the rain cloud. Previous 
methods were limited to fairly simple geometric shapes such as cubes and cylinders. 
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Abstract The last decade has seen tremendous growth in cloud dynamical and microphysical 
models that are able to simulate storms and storm systems with very high spatial resolution, 
typically of the order of a few kilometers. The fairly realistic distributions of cloud and 
hydrometeor properties that these models generate has in turn led to a renewed interest in the 
three-dimensional microwave radiative transfer modeling needed to understand the effect of 
cloud and rainfall inhomogeneities upon microwave observations. Monte Carlo methods, and 
particularly backwards Monte Carlo methods have shown themselves to be very desirable due to 
the quick convergence of the solutions. Unfortunately, backwards Monte Carlo methods are not 
well suited to treat polarized radiation. This study reviews the existing Monte Carlo methods 
and presents a new polarized Monte Carlo radiative transfer code. The code is based on a 
forward scheme but uses aliasing techniques to keep the computational requirements equivalent 
to the backwards solution. Radiative transfer computations have been performed using a 
microphysical-dynamical cloud model and the results are presented together with the algorithm 
description. 


1. Introduction 


Passive microwave remote sensing of rainfall has always been hampered by non-linear relations 
between rainfall rates and satellite observed radiances. Nearly 20 years ago, Weinman and 
Davies [1977] and then later Kummerow and Weinman [1988, a] built radiative transfer models 
to account for the 3-Dimensional nature of precipitation. Unfortunately, these early models were 
limited to simple cuboidal geometries and much of the 3-dimensional effect was determined by 
the assumptions regarding the size and shape of the clouds. More recently, cloud dynamical 
models such as Tao et ah [1987] or Tripoli [1992] have made great strides in generating fine 
scale cloud and storm properties which appear quite realistic when compared to ground based 
radars. Because the cloud models fully specify the atmospheric structure and hydrometeor 
contents of the storms, these models offer the next level of realism needed to understand the 
radiative effects of 3-dimensional cloud and rainfall fields. 

From a remote sensing point of view, the idea of using backward Monte Carlo simulations to 
deal with these new complex cloud structures in the microwave regime was due to Petty (1994), 
although that paper still dealt with abstract cloud structures. The work was quickly followed by 
Roberti et al., (1994), dealing explicitly with cloud dynamical model output and soon thereafter 
by Liu et al, (1996) who first considered polarized radiative transfer, although only in a limited 
sense. Backwards Monte Carlo methods were found to be very easy to implement, flexible and 
allowed an easy interpretation of the interactions of the radiation with the medium Large 


absorption cross-sections meant that photons were absorbed quickly in the microwave regime 
leading to numerical convergence with relatively few photons. With solution speeds comparable 
to analytical approaches, the Monte Carlo methods are indeed attractive due to their simplicity. 

One shortcoming of the backward Monte Carlo solutions is its inability to deal properly with 
polarization introduced by atmospheric constituents unless the vertical and horizontal extinction 
cross-sections are identical. Current satellite sensors such as the SSM/I and the TRMM 
Microwave Imager measure microwave brightness temperature (Tb) in both the horizontal as 
well as the vertical polarization. Over water backgrounds, emissivities are significantly different 
for these polarizations and the Tb differences may thus be used as a measure of the total 
attenuation of the microwave signal as it propagates through the cloud. Land surfaces generally 
produce only very small polarization signatures. These signatures are almost impossible to 
detect in the presence of rainfall which further attenuates any surface signal. An analysis by 
Heymsfield, and Fulton (1994), however, revealed 85 GHz polarization signals (TbH - TbV) of 
about 8-1 3K over the stratiform precipitation region of both extra-tropical and tropical 
Mesoscale Convective Systems (MCS) over land. There is a similar (even if not recurring) 
pattern at 37 GHz although the magnitude of the polarization signal does not appear to exceed 
7K. No polarization signals were found in the convective portions of these storms. Heymsfield 
and Fulton argue that these polarization signals cannot be attributed to the surface conditions but 
must be attributed to atmospheric constituents and speculate that they are probably due to non- 
spherical hydrometeors, such as smaller ice crystals, falling slowly with an approximate 
horizontal orientation in the less turbulent mesoscale updraft. The lower polarisation differences 



in the convective region could be caused by the presence of randomly oriented snow particles, 
due to the convective motion, and to the presence of bigger, almost spherical, graupel particles. 

Backward Monte Carlo schemes cannot, by design, deal with photons changing polarisation 
state. The main reason is that in a backward scheme the first collision is effectively the last 
collision in the temporal sequence. This implies that, when a photon is started, the initial 
polarization state is unknown. To overcome this problem, Liu et al., (1996) developed the so 
called "backward-forward" scheme. This scheme does allow for the treatment of polarized 
radiation, but requires that the extinction matrix be diagonal. Observational evidence such as the 
Heymsfield and Fulton study, which suggests preferentially oriented particles, cannot be dealt 
with using diagonal extinction matrices. The next development step for microwave radiative 
transfer codes must be fully polarized 3-dimensional methods. With this aim, a new forward 
Monte Carlo method has been developed for the solution of the radiative transfer equation to 
include a full treatment of the Stokes parameters. 

In the process of obtaining an efficient forward Monte Carlo code, four different Monte Carlo 
methods have been developed and are described in the following sections: 1) a 3-D backward- 
forward polarized code for spherical or randomly oriented hydrometeors similar to Liu et al., 
(1996); 2) a direct plane parallel non polarized code; 3) a direct plane parallel polarized code for 
spherical or random oriented spheroidal particles; and 4) a direct plane parallel polarized code 
for oriented spheroidal particles. The last section examines the necessary extensions to make a 
fully polarized version work within a 3-dimensional framework. The aim of these codes is to 
simulate the radiance that could be measured by a space-bome radiometer at the top of the 



atmosphere in a specific direction needed to build a quantitative understanding of the effect of 
non-sphericai particles upon the polarization signatures. The Goddard Cumulus Ensemble 
(GCE) model Tao et al. [1987] serves as a realistic hydrometeor field which can be used to study 
the behavior for varied hydrometeor shape and density assumptions. 

2. The cloud model and the computation of scattering parameters. 

All the radiative transfer computations in this study refer to realistic cloud profiles obtained form 
the cloud microphysical-dynamical model developed by Tao et al. [1987], The cloud domain 
consists in 64x64 pixels at 15minutes in the evolution of the simulated storm. A single time has 
been chosen which represents a mature squall line. The top boundary is considered as cosmic 
background with a temperature of Tc=2.7K, while the surface has been assumed to be 
Lambertian with a temperature of 300K and emissivity e s = s h = e v = 0.9 . Each pixel consists 
of 21 vertical layers from the surface to a height of 9 km. The cloud model specifies the height, 
pressure, temperature and relative humidity of each layer. Within each layer, cloud water, rain 
water, cloud ice, snow and graupel contents are specified. Graupel is assumed to be spherical 
and its density is 0.4 gm cm' 3 . Both rain and snow are assumed to be horizontally oriented 
oblate spheroids. The asymmetry factor for rain is the same used in Kummerow and Weinman 
[1988,b] and its density is 1 gm cm' 3 . A Marshall and Palmer size distribution is used for snow, 
rain and graupel. Tw r o different codes have been used for the computation of the scattering 
parameters. The first uses the T-matrix method {Barber and Hill [1990]) while the second uses 
the Discrete Dipole Approximation {Draine and Flatau [1988]). The codes by Barber and Hill 
[1990] allow the computation of the scattering matrix only for randomly oriented particles, but, 


for the considered cases, these codes were faster and more accurate of the DDA codes. It is 
reasonable to suppose that the polarization differences are mainly due to the different scattering 
and absorption coefficients and therefore, in the following simulations, the matrices for 
randomly oriented particles have been used. The obtained results could be slightly altered by this 
simplification, but the physical considerations are still valid. For spheroidal oriented particles, 
the extinction coefficient and the single scattering albedo have been computed both for the 
vertical and horizontal polarization for 30 different polar directions of incidence between 0° and 
90°. The particle symmetry allows no computations to be done for different azimuth directions. 


3. Polarized 3-D backward-forward Monte Carlo method for spherical or random 
oriented hydrometeors 

A backward 3-D Monte Carlo method previously developed by Roberti et al. [1994] has been 
extended to allow the computation of polarized radiation. Here only the necessary modifications 
are described. A similar Monte Carlo method has been developed by Liu et al. [1996], In 
particular, here and in the following, we adopt the notation of Chandraseckar [ 1 960] describing 
a beam of polarized radiation with the Stokes parameters I = (I h ,I v ,U,V) where h and v refer 
to two mutually orthogonal directions respectively parallel and perpendicular to the reference 
plane. In our case the reference plane is the XY plane. The photons are released at the point of 
each subcloud where the T B is to be computed, with the direction opposite to the one in which 
they would physically propagate. At each scattering event the position of the photon, as well as 



the incident and scattering directions are memorized. When the photon is finally absorbed, then 
the initial Stokes parameters are 1° = (T X T X 0,0) where T x is the physical temperature of the 
medium at the point of absorption. This assumption canbe made because the Planck function is 
virtually linear with temperature in the microwave regime for typical atmopheric temperatures, 
and can therefore be replaced by T x . This further implies that the radiances can be interpreted as 
brightness temperatures rather than power per unit area. Then the path of the photon is traced 
forward through the scattering events and the scattering matrix is computed at each step. The 
Stokes vector after a scattering event is obtained from the Stokes vector before the scattering 
(primed variables) with the transformation: 
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The definition of angles /, and i 2 is consistent with the one given by (Chandrasechkar [1960], 
pp.39) and the scattering matrix P coincides with the definition given by ( Van de Hulst [1957], 
pp.44). In particular P u is the phase function for the unpolarized radiation. The division by P n 
in Eq. (1) is necessary to remove the bias introduced by the biased sampling of direction after 
scattering. In fact rj = cosi9 is sampled from P u and r, is sampled randomly between 0 and 2n , 
which are only a first approximation of the correct distributions ( Kattawar and Plass [1968], 
Collins et al. [1972]). The angle i 2 can computed from trigonometry (Lion [1980] pp.223). The 
final Stokes vector is finally derived as 


1 = S n ...5,I 0 (5) 

where N is the number of scattering events. The surface can be either specular or Lambertian 
with a given angular and polarization dependent emissivity. If an interaction with the surface 
occurs, a random number r uniformly distributed between [0,1] is compared with both the 

vertical and horizontal emissivities, sfj and e v . Three different events may occur. 

If 


r <s h , r<e v 


( 6 ) 



then the photon is absorbed and the corresponding Stokes vector is I 0 =(T S ,T S ,0,0) , where T s 
is the surface temperature. If 


r>s h , r > £ v 


then the photon is scattered, else if 


(7) 


r<e h , r>s v (8) 

then the photon is scattered and traced through the successive scattering events until it is finally 
absorbed with I 0 = (0. T x ,0,0) . Then the scattering events are traced forward and the Stokes 
vector is multiplied by the appropriate scattering matrices as in Eq. (1), until the interaction with 
the surface occurs. At that point the surface temperature T s is added to the I h component of the 
current Stokes vector. If the surface is Lambertian the scattering matrix is given by 
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An analog procedure is followed if r > sfr , r< e v To test the validity of the code, the results 
have been compared with those in Weinman and Guetter [1977] for a plane parallel cloud and 
they have been found in good agreement. This code can successfully take into account spherical 


or randomly oriented spheroidal particles, in which cases the extinction and scattering 
parameters are equal for the two polarization. Only the surface emissivity can be different for the 

two polarization. The code can take into account horizontal as well as vertical inhomogeneities 
of the medium. 

4. Non polarized direct Monte Carlo method 

The inability of backward MC, or even backward-forward MC methods to deal with non- 
diagonal scattering matrices requires the development of a polarized forward MC method. The 
first step was to develop a non-polarized forward Monte Carlo code competitive in terms of 
speed of convergence with the backward Monte Carlo code developed by Roberti et al. [1994], 
For this and the following codes the medium is considered to be plane parallel, even if the 
method could be extended to consider the horizontal inhomogeneities of the cloud as suggested 
in Sec.6o. With a forward Monte Carlo method the photons are released by the medium 
coherently with the source parameters and their propagation is followed until the escape from the 
medium. Unfortunately, to simulate the T B that could be measured by a space-borne radiometer, 
only the photons which escape the medium in a specific point and direction will be taken into 
account. This usually involves high computational times because only a fraction of the photons 
contribute to the radiance. To overcome this problem a number of biasing techniques (see also 
Roberti [1997]) are used. 


4.1 Source parameters 


The first problem in the definition of the source concerns the choice between the atmosphere, the 
surface, or the cosmic background as an emission source and has been solved by computing the 
irradiance which is proportional to the number of photons emitted. The total irradiance for a 
plane parallel atmosphere of height Zg is equal to 


F a - \ fs(z)T(z)k(z)dQdz = 4n js(z)T(z)k(z)dz (10) 

0 4rc 0 

where e,T, k are respectively the emissivity, the physical temperature and the extinction 
coefficient at height z. The irradiance of a Lambertian surface is equal to 

F s = \ £ X COS0 dQ = ne s T s (1 1) 


where s s and T s represent respectively the surface emissivity and temperature. The irradiance of 
the cosmic background of temperature T c =2.7 K is equal to 


F c = it T c 


( 12 ) 


The photon is released from the atmosphere if 


( 13 ) 
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and in a similar manner for the other sources. In the atmosphere a uniform emission in the 
vertical direction has been chosen. The direction of emission will be isotropic for the atmosphere 
and Lambertian for the surface. 

4.2 Tracing procedure 

The distance to collision is computed using two different methods depending on whether the 
photon is downwelling or upwelling. If the photon extended path intersects the surface, the 
optical distance to collision r coU is computed as 

T con = “ In M (14) 


where r is a random number uniformly distributed between [0,1]. If the photon extended path 


intersects the upper boundary of the atmosphere, a collision is forced before the photon escapes 
the. medium by selecting the optical distance from the truncated exponential distribution Collins 
et al [1982] so that 


T con =-ln[l-r(l-<T r ')] 


(15) 


where r, is the optical distance along the photon path, from the current photon position to the 

upper boundary of the atmosphere. In this case the photon weight is adjusted multiplying it by 
the factor 

0-e" r ) (16) 

After computing the distance to collision the photon can find itself in two possible situations. 

1 ) If the photon crosses a boundary before interacting with the medium it is advanced to the 
boundary and a new distance to collision is computed taking into consideration the distance 
already traveled House and Avery [1969], The photon is then relaunched from the boundary it 
would have crossed. 

2) If a collision occurs, then a scattering event is forced, thus avoiding absorption. The 
corresponding bias is removed by multiplying the photon weight by the single scattering albedo 
A (z) of the medium in the current position of the photon, which represents the probability of 
photon survival. If the collision is with the surface the photon weight is multiplied by (1 — e s ) . 
The new direction after scattering is computed using the phase function. Furthermore, to reduce 
computational times, the contribution that the photon would given to the T B if scattered in the 
viewing direction and if it would propagate to the upper boundary of the atmosphere without any 
further interaction with the medium is computed and the photon weight is modified accordingly. 
In the following W jn denotes the weight effectively attached to the n-th photon during the 
propagation at the i -th scattering event and W’ in indicates the weight attached to the n-th 
photon when its contribution to the T B if scattered (or emitted) in the viewing direction is 
computed. The weights are modified as follows 
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where A,, is the single scattering albedo at the considered layer, /?(cos0,) is either the 
atmospheric or the surface phase function evaluated on the angle between the incidence and the 
viewing direction, r/ is the optical depth from the point in which the /-th collision has occurred 
to the upper boundary along the viewing direction. The phase function /?(cos0 o ) for the 
emission contribution is l/4n for the atmosphere and cos(^)/n for the Lambertian surface where 
9 V is the polar viewing angle. The photon continues to propagate until the associated weight 
becomes less than the threshold below which the photon is eliminated. The threshold value has 
been fixed at 10E-6 in all the computations that follow. 


4.3 T b computation 

The irradiance emitted by the atmosphere is computed as the sum of the contributions of all the 
layers of infinitesimal thickness dz (Eq.(10)). With a Monte Carlo method, the integral in 
Eq.(10)is computed numerically as 


(18) 
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where is the number of photons emitted by the atmosphere and Z 0 is the atmosphere height. 

s mT„,k n depend on the height of emission z of the n-th photon. The fraction F a of F a which 
effectively contributes to the T B is 
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where N is the number of scattering events of the considered photon (/= 0 indicates emission), 
p( cos©) e r/Q v represents the probability that the photon after a scattering event propagates to 
the upper boundary of the medium in the viewing direction (i.e. in an infinitesimal solid angle 
df> x ) without undergoing any other collision. The atmospheric contribution to the T B can be 
computed by dividing F n by cos(6y)dO v . 
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The contribution to the T B given by the surface emission is 


( 21 ) 
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where A(y is the number of photons emitted by the surface. The contribution to the T B given by 
the cosmic background is 
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where A/,- is the number of photons emitted by the cosmic background. Consequently the T B can 
be computed as a sum of Eqs. (20), (21) and (22). 


T, = + T„, + T tc (23) 

The results of this code have been compared with those of an equivalent backward code [Roberti 
et al, 1995]. In Fig. 1 the resulting T B are presented both for a subcloud in the anvil (stratiform 

region) and a subcloud in the core (convective region) for 3 different frequencies. The 
convergence of the forward code (dashed line) is slower and the oscillations are higher. 
Nevertheless, in the worst case, only 500000 photons are needed to keep the oscillations below 
IK. 


5. Direct polarized Monte Carlo method (spherical or random oriented spheroidal 
hydrometeors) 

The code described above has been extended to include the polarization to treat spherical or 
random oriented spheroidal hydrometeors. 

5.1 Source Parameters 

The choice of the photons emission source is based, as for the unpolarized radiation, on the 
irradiances of the different sources. For spherical or random oriented spheroidal particles the 
extinction matrix for the Stokes vector is diagonal K= diag(k) which means that the four 
components of the Stokes vector are subject to the same extinction. Therefore the distance to 
collision can be computed as for the unpolarized radiation. 

5.2 The tracing procedure 

After emission, a scattering event is forced at every collision by multiplying the photon weight 
by the single scattering albedo and the new direction is determined by sampling the phase 
function. The new Stokes vector after a collision is computed as in Eq. (1). As explained in Sec. 

4.2 a drastic reduction of computational times is obtained by considering the contribution to the 
final Stokes vector given by each photon collision. The photon weight of the n-th photon is 


modified as follows: 
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The matrices S in Eq.(24) are evaluated on the angle between the photon incidence and the 
scattering direction, while the matrices Sare evaluated between the incidence and the viewing 
direction. The expression for matrices S is given in Eq.(l) and 


t j- i * a A* + A A a + A 
L =diag\ A„,A V , , — 1 


(25) 


For spherical or random oriented aspherical hydrometeors the single scattering albedos for the 
horizontal and vertical polarization A a , A v are equal. The photon continues the propagation until 

the first two components of the Stokes vector become lower of a fixed threshold and the photon 
is eliminated. 


5.3 Computation of the Stokes vector 


As for the unpolarized radiation the contribution to the final Stokes vector given by photons 
emitted respectively by the atmosphere and the surface (cosmic background ) is respectively 
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The Stokes vector at any computational point and direction is given by 


I=I „+Is+Ic (28) 

Numerical results will be presented in Sec.6.4. 

6. Direct polarized Monte Carlo method for oriented spheroidal particles. 

Applying a backward or even a backward forward Monte Carlo scheme when oriented 
spheroidal particles are involved gives rise to a number of problems. The angular dependence of 
the single scattering albedos is such that the probability of scattering is different depending on 
the photon propagating towards or from the detector. The use of a direct method seems ideal to 
avoid such problems. Furthermore the extinction matrix is of the type 
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and therefore a different extinction coefficient is associated to each Stokes component and this 
must be carefully considered in the photon propagation. Since the matrix is not diagonal, the 
extinction of each Stokes component depend also on the values of other components. A 
backward method seems to be useless in this respect because the Stokes vector is unknown 
during the propagation. 


6.1 Source parameters 


The choice of the emission source is made in a way similar to the ones previously described. 
For the atmosphere the following quantities are computed, which represent the average number 
of photons emitted by the atmosphere and the surface: 
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F a h > F a v > F s h ’ F s v are Oie same quantities as in Eqs. (10-11) except that the emissivities can be 
different for the horizontal and vertical polarization. Eq.(12) is still valid for the cosmic 


background. These quantities are necessary for the computation of the unknown Stokes vector. 
To increase the speed of convergence, it is useful to employ quantities which represent the real 
proportion between the energy emitted by the different sources. In the following the extinction 
matrix will be considered as diagonal, neglecting the circular polarization term. It is in fact 
possible to assume, in the case of thermal emission, that U and V are relatively small with 
respect to /*and /,. To overcome the problem of different extinction coefficients for the 
different components of the Stokes vector, another biasing techniques is introduced sampling the 
distance to collision from a biased probability distribution and multiplying the Stokes vector by 
proper coefficients to remove the statistical perturbation introduced. For downwelling photons 
the biased probability density for the optical distance to collision is 

p(z)dr = e~ kx d(kx) (32) 


where k is an extinction coefficient computed as 
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where kf, and k v are the extinction coefficients for the two polarizations. This choice has been 
made in order to minimize the multiplying factors for the elimination of the bias. The distance to 
collision is given by 
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and the matrix for the correction of the bias is given by 
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For upwelling photons, as explained in Sec. 4.2, a collision is forced before photons escape the 
medium and the distance to collision becomes 


d coii =-ilog[l-r(l-e~ ? )] (36) 

where f is the optical thickness computed from the current photon position to the top of the 

atmosphere in the direction of propagation using the extinction coefficient £. The weight 
correction matrix is similar to the one in Eq.(35), except that each element is still multiplied by 
a factor (1 — e r ). Some attention must be paid to the computation of the new distance to 
collision after the crossing of an internal boundary of the medium. The true distance to collision 
probability densities for a photon crossing the boundary between layer 1 and 2 are 


P h (r)dT = Q~ k,l ' d 'Q k "^ d - d ^d{k hl d) d>d ] 
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for the H and V polarized photons, where d x is the distance between the current photon position 
to the boundary in the direction of propagation of the photon. The biased probability density 
with which the distance to collision inside the new layer is corrected taking into account the 
distance already travelled is 
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The distance to collision in the second medium becomes 
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where d coll _ olJ is the distance to collision computed in the first medium. Therefore the weight 
correction matrix is a diagonal matrix with elements 
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C is equal to (1 - e r ) for upwelling photons and C =1 in the opposite case. 


6.2 The tracing procedure 


The photon weights during the propagation are modified as follows 
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^V.’=S (>1 K -\z M e t ) K(z i+i 0J E, +I W, 


where 



E ( . = diag(t' T '" ,e‘ r '-' ,e' rui ,e’ rv ‘ ) 


(42) 


a nd T(j hV jj,V),i are the optical thickness from the /-th scattering point to the observation point 
computed with the appropriate extinction coefficient 


K (. z,9 ) = cling k h (z,9),k v (z,9 ) 


k h {z,9) + k v (z,9) k h {z,9) + k v {z,9) 


(43) 


The multiplication by K' 1 in Eqs. (41) is introduced to take into account the angular dependency 
of the extinction coefficient. Matrix L is Eq.(25). The final Stokes vector is computed with 
formulae analogous to Eqs. (26-28). 


6.3 Extension of the code to the 3-D case 


The main problem in the extension of the forward codes to the 3-D case deals with the treatment 
of the emission sources. If the horizontal inhomogeneities of the medium are to be taken into 
account photons should be released by the different layers of the different subclouds according 
to the thermal contribution. For each subcloud nx, ny a threshold must be fixed as (see Eqs. 
10 , 11 , 12 ) 


Trs 


nx,ny 


+ F +F 


(44) 



If NXCL, NYCL are the subcloud numbers in the x and y directions, a photon is emitted by 
subcloud hx, ny if a random number r uniformly distributed between [0,1] satisfies 


Z Trs >.i 

i=\,rtx j = \,ny - 1 

Z Trs ,i 

/ = 1 , NXCL j=\,NYCL 




Z Trs ,j 

i=\,nx j = \,ny 

Z Trs u 

i=\,NXCL j=\,NYCL 


(45) 


Once the subcloud has been determined, the vertical position for emission inside each subcloud 
can be prescribed as in Sec. 4.1. If the subcloud horizontal dimension is small compared to the 
vertical dimension, the photons are simply started at the center or the subcloud. Emission 
points must be sampled uniformly in the case that the subcloud dimension is large compared to 
the vertical dimension. The computation of the Stokes vector must be partially modified since 
the Stokes vector must be computed at the top of each subcloud. At each scattering event a line 
must be ideally drawn from the point of collision to the top of the cloud in the viewing direction 
The line will intersect the top of only one of the subclouds. The scattering event will give a 
contribution only to the Stokes vector associated with that subcloud and the contribution can be 
computed as described in Sec. 4.3 for the unpolarized case. The rest of the algorithm remains 
unchanged. Photons crossing the vertical boundary between two subclouds or the horizontal 
boundary between two layers are treated analogously to the plane parallel methods. 


6.4 Numerical Results 


The convergence of the two forward polarized codes has been analyzed with reference to the 
backward-forward code. The resulting J H and A T V _ H are presented for a subcloud in the 
stratiform region (Fig. 2) for 3 different frequencies. Similar results can be obtained for a 
subcloud in the convective region. A Lambertian surface with emissivities e t = e h =e v = 0.9 
has been chosen in order to analyze the polarization signature introduced by atmospheric 
scattering. Considering spherical hydrometeors the resulting A T y _ H are negligible. The same 
considerations done for the forward unpolarized code are still valid here. The convergence of the 
forward codes is slower, but, if the required accuracy is of the order of IK, an average number of 
500000 photons is sufficient. 

To simulate the observed polarized results ( Heymsfield and Fulton [1994]), a cross section of 
the cloud model, made of 64 subclouds, normal to the main convective line is considered. A 
sensitivity analysis of the A7^,_ H as a function of snow density, asymmetry factor (a/b) and 
concentration using plane parallel computations on single subclouds is performed. In order to 
avoid changing the total ice concentration, snow concentration is varied by reducing graupel and 
keeping the total ice mass constant. The results are presented for a subcloud in the stratiform 
region (Fig. 3) and one in the convective region (Fig. 4). The resulting AT y _ H are shown, as a 
function of the asymmetry factor, for density of snow equal to 0.1 and 0.2 g/m\ 5 different 
percentages of graupel transformed to snow (Perc=0, 0.25, 0.50, 0.75 and 1.00) and at 3 
frequencies (19, 37 and 85 GHz). In particular the concentrations of snow and graupel have been 
computed as 


snow = snow + graupel * Perc 


(46) 


graupel = graupel * (1-Perc) 

The general behavior of the plots is the one expected: the A T V _ H increase at increasing density 
and percentage of graupel transformed to snow. Furthermore the AT y _„ are higher when the 
asymmetry factor is higher, and this appears more clearly when the percentage of graupel 
transformed to snow is more than 50%. If the snow content is low (Perc=0, 0.25), the A T v _„ at 
lower frequencies are higher than those at 85 GHz. This might be due to the fact that the 
radiation at 85GHz is more sensitive to the depolarizing effect of graupel scattering. To obtain 
the desired A T V _ H in the anvil the content of snow originally present in the cloud seems to be 
too low. It is necessary to transform 50% of the graupel to snow to start obtaining a AT y _ H of 
1 IK at 85GHz. . On the other hand, if all the graupel is converted to snow, the A T V _ H become 
too high. In the core of the cloud the increase of snow content does not give AT V _ H similar to 
the observed ones. For Perc< 0.5 the influence of the remaining graupel is still too high, 
maintaining the resulting AT V _ H relatively low. On the other hand, for high content of oriented 
spheroidal snow particles [Perc> 0.75), the AT y _ H become too high at all frequencies, especially 
for high values of the asymmetry factor. These considerations seems to confirm the hypothesis 
of the presence of bigger spheroidal or randomly oriented ice particles in the core of storm, due 
to the updraft motion. 


Since the percentage of graupel that, transformed to snow, gives A T v _„ in accordance with the 
observed results varies with the different subclouds, we tried to transform graupel to snow, 
starting from the higher layers of each of the 64 subclouds, for a total of 1.5kg/m 2 . The resulting 
r„and A T r _ ff are shown at 19,37 and 85 GHz in Fig.5A, B respectively. The density of snow is 
0.1g/cm\ the asymmetry parameter is 0.3. 

7. Conclusions 

Four different Monte Carlo codes for the computation of microwave radiation emerging from 
cloud structures have been presented. Three codes allow to take into account polarization. One 
of them is 3-D, but treats only spherical or randomly oriented hydrometeors which give only 
small polarization signal. The code which can deal with oriented spheroidal particles is currently 
plane parallel. Future work will concentrate on the extension of this code to the 3-D case. The 
convergence and accuracy of the forward Monte Carlo codes is lower with respect to backward 
codes, but computational times are limited thanks to the biasing techniques that have been used. 
Nevertheless Forward Monte Carlo codes they seem to be the only way to treat polarization 
when oriented aspherical particles are involved. 


One of the forward Monte Carlo codes has been used for a sensitivity analysis of the A T V _ H as a 
function of snow parameters. To summarize the results it is possible to say that, while many 
combinations of density and asymmetry factor can explain the observed AT y _ H , it was found 
that converting a total of 1.5 Kg/m 2 of graupel into snow throughout the cloud model 
simulation produced generally adequate agreement with observations. The results represent a 
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contribution to the understanding of clouds microphysics in relation to radiometer observables 
at millimeter frequencies. While the poor knowledge of the real cloud parameters is underlined, 
the effort is toward creating a foundation for realistic cloud modeling. 
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FIGURE CAPTIONS 


Figure 1 : Tg for a subcloud in the anvil region and one in the core region of the cloud for 
density of snow 0.1 g/cm 3 , 3 different frequencies (19.35, 37.0, 85.6 GHz.), obtained with the 
non polarised backward Monte Carlo code ( )and the forward code ( ) . 

Figure 2: Tfj and A Ty.}{ for a subcloud in the stratiform region of the cloud for density of snow 
0.1 g/cm 3 , 3 different frequencies (19.35, 37.0, 85.6 GHz.), obtained with the polarized 

backward- forward Monte Carlo code ( ), the forward code for spherical particles ( ) 

and the forward code for oriented spheroidal particles (• • •) . 

Figure 3: tsTy.jjiox a subcloud in the stratiform region of the cloud for density of snow 0.1 and 

0.2 g/cm 3 and 5 percentages of graupel transformed to snow (0, 25,50,75 and 100%).( ) 19.35 

GHz, ( ) 37.0 GHz, (• • ■) 85.6 GHz. 


Figure 4: A7y_//for a subcloud in the convective region of the cloud for density of snow 0.1 
and 0.2 g/cm 3 and 5 percentages of graupel transformed to snow (0, 25,50,75 and 100%). 
( ) 1 9.35 GHz, ( ) 37.0 GHz, (• • •) 85.6 GHz. 

Figure 5: 7// and A Ty.jj associated to the considered cloud section. The density of snow is 
0.1 g/cm 3 , the asimmetry parameter a/b is 0.3. ( )19.35 GHz, ( )37.0 GHz ,(• • •) 85.6 


GHz 












































